########################################################################
# PLS 206 — Week 5 live coding, part 1: scoring models
# Grey Monroe, Fall 2026
########################################################################

# install.packages(c("ggplot2", "MuMIn"))
library(ggplot2)
library(MuMIn)

# 50 farms, cover-crop biomass (kg/ha) and 15 candidate predictors
farms <- read.csv("https://greymonroe.github.io/PLS_206/data/cover_crop.csv")
farms <- subset(farms, select = -c(farm, clay))   # sand + silt + clay = 100
names(farms)

########################################################################
# 1. Overfitting: 25 farms vs the other 25
########################################################################

set.seed(5)
tr <- sample(50, 25)            # fit on these 25 farms
fit_on <- farms[tr, ]
new    <- farms[-tr, ]          # predict these 25

m1  <- lm(biomass ~ gdd, data = fit_on)
m15 <- lm(biomass ~ ., data = fit_on)

summary(m1)$r.squared           # 0.57
summary(m15)$r.squared          # 0.88: R-squared only goes up

rmse <- function(obs, pred) sqrt(mean((obs - pred)^2))
rmse(new$biomass, predict(m1, new))    # 1164 kg/ha
rmse(new$biomass, predict(m15, new))   # 2573 kg/ha: much worse on new farms

par(mfrow = c(1, 2))
plot(new$biomass, predict(m1, new), main = "gdd only", xlab = "Observed", ylab = "Predicted")
abline(0, 1, lty = 2)
plot(new$biomass, predict(m15, new), main = "all 15", xlab = "Observed", ylab = "Predicted")
abline(0, 1, lty = 2)
par(mfrow = c(1, 1))

########################################################################
# 2. Likelihood and AIC by hand
########################################################################

m_temp <- lm(biomass ~ gdd, data = farms)

# each farm: how high it sits on the model's bell, as a log
sigma <- sqrt(mean(resid(m_temp)^2))                   # width of the bell
y <- farms$biomass[4]; y_hat <- fitted(m_temp)[[4]]    # one farm
1 / (sigma * sqrt(2 * pi)) * exp(-(y - y_hat)^2 / (2 * sigma^2))
dnorm(y, mean = y_hat, sd = sigma)                     # the same
scores <- dnorm(farms$biomass, mean = fitted(m_temp), sd = sigma, log = TRUE)
round(range(scores), 2)         # worst and best farm
sum(scores)                     # add them up
logLik(m_temp)                  # -420.47, df = 3: intercept, slope, sigma

-2 * as.numeric(logLik(m_temp)) + 2 * 3
AIC(m_temp)                     # 846.93: the same

# the model, against the data
ci <- data.frame(gdd = seq(min(farms$gdd), max(farms$gdd), length.out = 100))
ci <- cbind(ci, predict(m_temp, ci, interval = "confidence"))
p <- ggplot(farms, aes(gdd, biomass)) +
  geom_point() +
  theme_classic()
p
p + geom_ribbon(data = ci, aes(gdd, ymin = lwr, ymax = upr), inherit.aes = FALSE, alpha = 0.3) +
  geom_line(data = ci, aes(gdd, fit), colour = "blue")

# the null model: the mean and sigma, k = 2
m_null <- lm(biomass ~ 1, data = farms)
logLik(m_null)
AIC(m_null)                     # 894.1: 47 worse

########################################################################
# 3. AICc, BIC
########################################################################

AICc(m_temp)                    # 847.46: + 2k(k+1)/(n-k-1)
BIC(m_temp)                     # 852.67: log(n) per parameter
log(50)                         # 3.9 per parameter instead of 2

########################################################################
# 4. Six ideas, model.sel()
########################################################################

m_null  <- lm(biomass ~ 1, data = farms)
m_temp  <- lm(biomass ~ gdd, data = farms)
m_moist <- lm(biomass ~ rain_total, data = farms)
m_soil  <- lm(biomass ~ organic_matter + soil_ph, data = farms)
m_tm    <- lm(biomass ~ gdd + rain_total, data = farms)
m_all   <- lm(biomass ~ gdd + rain_total + organic_matter + soil_ph, data = farms)

ms <- model.sel(m_null, m_temp, m_moist, m_soil, m_tm, m_all)
ms
round(as.data.frame(ms)[, c("df", "logLik", "AICc", "delta", "weight")], 2)

# nested: the F test agrees
anova(m_temp, m_tm)             # p = 0.62, delta AICc 2.11

# not nested: no F test, AIC still works
anova(m_moist, m_soil)

########################################################################
# 5. The two traps
########################################################################

# same response: log(y) and y cannot be compared
m_log <- lm(log(biomass) ~ gdd, data = farms)
AIC(m_temp, m_log)              # 847 vs 123: different units

# same rows: missing values change n
farms_na <- farms
farms_na$soil_ph[c(3, 17, 29)] <- NA
m_a <- lm(biomass ~ gdd, data = farms_na)
m_b <- lm(biomass ~ gdd + soil_ph, data = farms_na)
c(nobs(m_a), nobs(m_b))         # 50 vs 47
AIC(m_a, m_b)                   # m_b looks better only because n is smaller

# fix: drop incomplete rows once, then fit everything
farms_cc <- na.omit(farms_na)
AIC(lm(biomass ~ gdd, data = farms_cc), lm(biomass ~ gdd + soil_ph, data = farms_cc))